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Abstract 

In a series of papers by Smirnova and Hayakawa, fractal 
analyses of the ULF emissions in the frequency range of f = 
0.001-0.3 Hz have been performed based on the geomagnetic 
data obtained in seismic active regions including Guam 
Island and Izu Peninsula. A certain dynamics of the 
spectrum slopes and the corresponding first-order fractal 
dimensions have been revealed in relation to the preparation 
phase of some strong earthquakes (EQs). In order to advance 
such a fractal approach which seems to be very promising 
for development of the EQ forecasting methods, we consider 
multifractal aspects in the analysis of geophysical data 
dynamics. First, as a case study, the seismicity distribution in 
the Kobe area of Japan is concerned in relation to the 
powerful Kobe EQ on 17 January 1995. Applying a 
multifractal approach, we have revealed that there is a 
gradual decrease in the certain high-order fractal dimensions 
which were obtained from the spectrum of singularities 
(multifractal spectrum) of the local seismicity distribution 
when the date of the Kobe EQ was approached. Many 
examples of the multifractal spectrum dynamics in relation 
to major EQs is introduced as a statistical study. It is then 
concluded that the multifractal analysis of geophysical data 
could be a promising tool for extraction of the precursory 
signatures of the extreme natural events including strong 
EQs. 
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Introduction Motivation for a Multifractal 
Approach 

The concept of self-organized criticality (SOC), which 
has been first introduced by Bak et al. (1988) to explain 
the 1/f noise (flicker noise) and scale-invariant (fractal) 
structure, can be considered as a general principle 
governing the behavior of a certain class of complex 
dissipative dynamic systems toward a rupture. This 


SOC theory widely used to investigate the dynamics 
of natural hazard systems including earthquakes (EQs) 
(e.g., Goltz, 1997) has been utilized in this paper as one 
of the principal points in our complex approach to 
study the EQ preparation processes. 

A lot of evidence has been accumulated on the 
important role of electromagnetic precursors of EQs 
for the purpose of short-term EQ prediction (e.g., 
recent books by Hayakawa and Molchanov (2002), 
Molchanov and Hayakawa (2008), Hayakawa (2009, 
2012)). Since the first attempt of fractal analysis for the 
seismogenic ULF emissions (Hayakawa et al., 1999), 
a phenomenological model of the large-scale 
evolutionary processes between two violent EQs based 
on the SOC concept has been suggested. Four 
principal stages in the EQ preparation processes have 
been taken into consideration: initial phase (random 
chaos just after the first EQ when the tectonic energy is 
fully released), subcritical, critical and super-critical 
stages. Supercritical stage is the final stage of the SOC 
evolution, where there is a high probability of the next 
violent EQ. Subcritical and critical stages are 
intermediate stages of the SOC evolution. Since the 
principal feature of the SOC dynamics is a power-law 
distribution (or fractal organization) of the system 
parameters both in space (scale-invariant structure) 
and time (temporal-invariant structure), fractal 
methods can be used to investigate nonlinear scaling 
characteristics of such a distribution at different stages 
of the strong EQ preparation. In our previous studies 
(Smirnova, 1999; Gotoh et al., 2003, 2004; Smirnova et 
al., 1999, 2001, 2004, 2010; Ida et al., 2005, 2006; Ida and 
Hayakawa, 2006; Smirnova and Hayakawa, 2007; 
Hayakawa and Ida, 2008), the focus was on ULF 
(f=0.001-lHz) geomagnetic measurements in seismo- 
active regions, since those ULF waves are the most 
sensitive to the variation of geoelectric parameters of 
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Earth's crust at typical depth of EQs (see e.g., 
Hayakawa et al., 2007; Hayakawa, 2012). Then, it is 
concluded that precursory signatures of the strong 
EQs are manifested in dynamics of the scaling 
characteristics (spectrum slopes and fractal 
dimensions) of the ULF emissions registered in seismic 
active regions. Here the multifractal approach was 
introduced to extract the EQ precursory signatures 
from spatial and temporal variations of geophysical 
data. The primary question which arises, is "Why 
multifractal approach should be used?" The answer 
relates mainly to the clustering of seismicity and an 
analogy between the destruction processes and the 
process of a multifractal measure generation. In this 
paper, we do not deal with seismogenic ULF 
emissions, but with the simple EQ distribution. First of 
all, the seismicity forms spot-like structures in a wide 
range of scales as seen in Fig. 1 taken from 
Kiyashchenko et al. (2004). Such structures are similar 
to multifractal measure distribution. Second, the 
destruction process leads to the formation of 
multifractal-like structures. The simple mechanical 
analogy is illustrated in Fig. 2(a). When there appears a 
crack, the initial stress M is redistributed around its tips 
(Step 1). In the zone B, the stress level decreases. Let us 
assume that its average value is Mi. In the zone A, the 
stress increases to the value M 2 > Mi. On the second step, 
if the other cracks appear in the zones A and B, the stress 
field is redistributed in a similar proportion M 2 : Mi (see 
Step 2 in Fig. 2(a)). This process is very similar to 
the multiplicative process of multifractal measure 
generation (see the corresponding figures in Feder 
(1988) and Mandelbrot (1989)). 

CLUSTERING IN JAPANESE SEISMICITY 


1000 km 



FIG. 1 AN EXAMPLE OF SEISMICITY CLUSTERING ON 
DIFFERENT SPATIAL SCALES (SEE KIYASHCHENKO ET AL. 
(2004) FOR DETAILS). 

Multifractal Formalism 

Multifractal approach is a synthesis of the fractal 
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analysis, which is more general than the monofractal 
approach and can be presented by a spectrum of 
generalized fractal dimensions d(q) and spectrum of 
singularities f(a) (multifractal spectrum) (e.g., Goltz 
(1997), Badii and Roliti (1997), Hayakawa and Ida 
(2008)). The physical and mathematical meanings of 
these parameters will be given later. An example of 
such spectra for a binomial multiplicative process 
presented in Fig. 2(b) is shown in Fig. 3. Curves 1 
correspond to the case Mi=0.6 and M2=0.4, while 
Curves 2 correspond to Mi=0.8 and M2=0.2 (see 
Kiyashchenko et al. (2003) for details). It is seen from 
Fig. 3 that the range of generalized fractal dimensions 
d(q) is wider for the more contrast redistribution of the 
mass Mi, M 2 (more inhomogeneous media). The same 
feature can also be noted on the basis of the 
multifractal spectral behaviour (see f(a) curves in the 
right part of Fig. 3). The f(a) curve is known to extend 
with an increase in contrast of Mi and M 2 . A general 
case of the real multifractal spectrum is presented in 
Fig. 4. A lot of parameters taken from the f(a) curve 
can be informative: the Lipshitz-Holder exponents ao. 
Oil, Oimm, OCmax, the range Aa=amax-amin, the dimensions 
Do, d(l) = a\, etc. Therefore, by tracing f(a) changes as 
well as the changes in spectrum of the generalized 
fractal dimensions d(q), it is possible for us to study 
the dynamics of natural tectonic systems and to 
investigate EQ preparation processes. 

Multifractal Properties of EQs: World-wide 
Studies 

Up to now, many workers have studied fractal and 
chaotic features of EQs. So the hierarchical properties 
of seismicity have been firstly described in the early 
works of Sadovskiy (e.g. Sadovskiy et al., 1984). The 
more complete studies in this field are summarized in 
the Lecture Note by Goltz (see Goltz, 1997 and 
references therein). According to the paper by Haikun 
(1993) quoted in Goltz (1997) some precursory changes 
of the multifractal f(a) curve have been observed 
inrelation to the strong EQs occurring in the Datong 
area, China during the 1970-1991. The schematic 
illustration of such precursory behaviour is shown in 
Fig. 5. The following peculiarities observed before 
major EQs can be noted based on this illustration. 

• The multifractal spectrum is found to extend 
(corresponding A a increases). The increase in 
Aa means a transition from homogeneous 
(random, space filling) to heterogeneous 
(ordered, complex, clustered) patterns. 

• amin and amax shift to the right and ao increases. 
The shift to the right of amin means that the 
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clustering within the most clustered areas 
becomes more intense (the local fractal 
dimension increases within these vicinities). 
The analogous shift of amax indicates that the 
clustering within the sparse areas becomes also 
enhanced. Thus an overall increase in the 
degree of clustering can be expected. 

• The increase of ao shows that most clusters 
now possess a higher local fractal dimension 
than before. An increase of f(ao) means a rise of 
Do-capacity dimension which is not usually 
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observed. 

• The increase in Af(a) indicates that a change 
has occurred in the ratio of highly clustered 
and sparsely populated areas. 

There is another paper by Hirabayashi et al. ( 1992 ), in 
which the dynamics of d(q) is presented in Fig. 6 for 
the California (left part) and Japan (right part) EQs. 
For each (a) or (b), the left panel in the bottom refers to 
the seismic quiet condition and the right, the seismic 
active condition. 
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FIG. 2 THE QUALITATIVE ANALOGY BETWEEN THE DESTRUCTION PROCESSES (a) AND THE PROCESS OF A MULTIFRACTAL 

MEASURE GENERATION (b) 



FIG. 3 SPECTRUM OF THE GENERALIZED FRACTAL DIMENSIONS (a) AND SPECTRUM OF SINGULARITIES (MULTIFRACTAL 
SPECTRUM) (b) FOR A BINOMIAL MULTIPLICATIVE PROCESS PRESENTED IN FIG. 2(b). (CURVE 1: Mi=0.6, M 2 =0.4; CURVE 2: Mi=0.8, 

M 2 =0.2). 




FIG. 5 PRECURSORY BEHAVIOUR OF MULTIFRACTAL 
SPECTRA AS OBSERVED BY HAIKUN (1993). THE /fa) - a CURVE 
CHANGES FROM THE SOLID ONE TO THE DASHED ONE 
SEVERAL YEARS BEFORE A MAJOR EQ OCCURS 
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FIG. 6 CHANGES OF THE GENERALIZED FRACTAL DIMENSIONS d(q) FOR CALIFORNIA (a) AND JAPAN (b) EQS. FOR EACH (a) OR (b), 
THE LEFT PANEL IN THE BOTTOM REFERS TO THE SEISMIC QUIET CONDITION, WHILE THE RIGHT, SEISMIC ACTIVE CONDITION 

(HIRABAYASHI ET AL.,1992). 


Our Studies of Multiracial EQ Precursors 

The world-wide studies of the multifractal properties 
of seismicity distribution seem to be not so rigorous 
and complete, and their results are often inconsistent 
with each other. So we have carried out our own 
researches to understand if there are any pronounced 
precursory signatures in dynamics of the seismicity 
multiracial characteristics. 

Data Used for Multiracial Analysis of the 
Seismicity 

Two EQ hypocenter catalogs have been utilized for 
our multifractal analysis: (1) Japan University 

Network EQ Catalog Hypocenter s, published by 
Earthquake Research Institute, University of Tokyo 
(http://www.eri.u-tokyo.ac.jp). This catalog contains 
the EQs with magnitude M > 2.0 occurring in the area 
with geographical coordinates cp = 26- 48° N and A = 
128 - 148°E during the temporal period from 1985 to 
1996. (2) Southern California EQ hypocenter catalog 
containing the EQs with magnitude M > 1.0 occurred 
in the area with geographical coordinates cp= 1 - 60°N 
and A = 26 - 178°W during the temporal period from 
1982 to 1999 (http://www.scecdc.scec.org/catalog). 

Our consideration was restricted only to those events, 
which occurred at the depth less than 60 km, and the 
strong EQs were selected with magnitude M > 6.5 
from both catalogs to study the dynamics of regional 


seismicity distribution prior to those strong events. For 
each of the strong EQs considered, the sub-catalog of 
the seismicity registered in its surrounding area was 
selected for analysis. The surrounding area is taken in 
the form of a box with the side A (Fig. 7) centered at 
the EQ epicenter. The size A of the area has to be 
selected according to the size of the rupture area of the 
strong EQ. Also this area should contain sufficient 
amount of data for multifractal analysis of the 
seismicity in a sliding temporal window. The available 
size A of the surrounding area has been roughly 
estimated. The size L of the rupture area can be 
estimated using the empirical relation (Sobolev and 
Ponomarev, 1999; Sobolev, 1990): IgL = 0.6M - 2.5, 
where M is the EQ magnitude. Computational results 
of the L values are shown in Table 1. 



earthquake 

epicenter 


FIG. 7 BOX - LIKE AREA NEAR AN EQ EPICENTER. 

The effective radius of the surrounding area (equal to 
A/2), which has to be selected for analysis, must not 
exceed several lengths of the rupture of an expected 
strong EQ. 
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TABLE 2 THE PARAMETERS OF THE STRONG EQS SELECTED FOR ANALYSIS. 


N 

data 

<P 

\ 

A 

Z- km 

M 

S , knr 

1 

27/10/88 

36.3° N 

141.749° E 

40 

7.7 

400 x 400 

2 

02/1 1/89 

39.736° N 

143.39° E 

0 

6.6 

400 x 400 

3 

07/02/93 

37.658° N 

137.309° E 

26.5 

6.6 

400 x 400 

4 

04/10/94 

43.383° N 

147.928° E 

46.2 

8.1 

400 x 400 

5 

28/12/94 

40.434° N 

143.867° E 

14.1 

7.5 

400 x 400 

6 

07/01/95 

40.232° N 

142.425° E 

55.9 

7.2 

400 x 400 

7 

17/01/95 

34.583° N 

135.027° E 

33.4 

7.2 

400 x 400 

8 

1 7/02/96 

37.306° N 

142.643° E 

59.1 

6.7 

400 x 400 

9 

19/10/96 

31.798° N 

131.972° E 

44.6 

6.6 

400 x 400 

10 

24/1 1/87 

33.01° N 

1 15.85° W 

1 1.2 

6.6 

200 x 200 

1 1 

28/06/92 

34.2° N 

1 16.44° W 

1 

7.3 

200 x 200 

12 

16/10/99 

34.594° N 

1 16.271° W 

0 

7.1 

200 x 200 


TABLE 1 THE RUPTURE SIZE L OF AN EQ WITH MAGNITUDE M 


M 

4.0 

4.5 

5.0 5.5 

6.0 

6.5 

7.0 

7.5 

S.0 

LJin 

0.8 

1.6 

3.16 6.3 

12.3 

25.1 

50.1 

100 

200 


That is an important condition for searching the 
precursory phenomena using the methods suggested 
in this paper. In fact, if the size of the rupture area of 
the selected strong EQ is much smaller than that of 
surrounding area (for example, if the size of 
surrounding area is A = 400 km and the size of rupture 
area is L = 6.3 km for the EQ with M = 5.5), the main 
part of the seismicity located in this area has no 
influence on the process of preparation of that strong 
EQ. Since we have selected the EQs with M>6.5, the 
corresponding rupture area is L>25 km (see Table 1). 
In such a case, the size of the surrounding area A= 
200-400 km seems to be suitable for multifractal 
analysis: on the one hand, it is in appropriate relation 
with the rupture size; on the other hand, such area 
contains sufficient amount of data for the analysis of 
multifractal characteristics of seismicity distribution in 
sliding temporal windows of reasonable length. 
Taking into account the amount of data in the sub- 
catalogs used for analysis and the chosen temporal 
window (less than approximately 3 years in our case), 
the size of surrounding area 400 * 400 km 2 seems to be 
acceptable for Japanese EQs, and the size 200 x 200 
km 2 is suitable for the EQs of Southern California. We 
can take the smaller size of surrounding area for the 
EQs of Southern California since the Californian 
catalogs contain a larger amount of data. Overall, 12 
strong EQs with M > 6.5 have been selected from both 
catalogs. The information concerning these EQs is 


summarized in Table 2. The corresponding sub- 
catalogs of seismicity have been analyzed in the EQ 
surrounding areas of the size S km 2 . Here N is the 
serial number of EQ, cp, latitude of the epicenter. A, its 
longitude, z, the depth of the hypocenter, and S is the 
square of surrounding area centred on the epicenter of 
the EQ. 

Analysis Procedure 

For each of the sub-catalogs, the following multifractal 
characteristics have been calculated as seen in Fig. 4: S 
=d( 1) (information dimension), d( 2) (correlation 
dimension), amin, =d(+oo), and aox (the abscess of the top 
of the multifractal spectrum). 

Calculations have been fulfilled with a sliding temporal 
window, whose length was selected individually in each 
case to contain an adequate set of statistics (about 1000 
events). In our case, 2D (two dimensional) distribution 
of the seismicity hypocenters has been studied as seen 
in Fig. 8. The distribution function pi (A) is defined as: 



FIG. 8 EQ EPICENTERS CONTINUUM. 

N, 

p^=if a) 

where Nj is the number of EQs occurring in the cell 
with number j and size A, and N is the total number of 
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EQs in the considered region. The generalized fractal 
dimension d(q) of the order q can be determined from 
the Renyi entropy (Kiyashchenko et al., 2003). 

Let us introduce the Renyi entropy of order q: 


In V jP ] A 

In (/ (A)) = , if q*l 

1 -q 

In (/ 9 (A)) = -X^(A)ln( A (A)), if q=l. 


Then, the generalized fractal dimension of the order q 
is defined as: 


d(q) = lim 

A— >0 


ln(/,(A)) 
ln(l / A) 


(4) 


So the generalized fractal dimension d(q) can be 
determined as a slope of the best fit line representing 
the Renyi entropy I(q) versus 1/A in the log-log plot. It 
is possible to show that the Holder exponent a and the 
/ (a) curve can be retrieved from the spectrum of 
generalized fractal dimensions d(q) using the 
transforms: 


dq ^ 

f(a(q)) = q-d(q)(q- 1) 


Dynamics of four multifractal characteristics, namely, 
the minimal value of the Holder exponent OCmin — d(+ ao), 
the information dimension S = ai = d(l), the correlation 
dimension d( 2), and the abscissa of the top of 
multifractal spectrum aox have been studied prior to 
strong EQs listed in Table 2. These multifractal 
characteristics provide us with the important 
information on the inhomogeneity of seismicity 
distribution and the level of seismicity clustering in a 
wide range of scales. So, the entropy S can be treated as 
the measure of inhomogeneity of the distribution of the 
seismicity: the lower (higher) values of S correspond to 
more (less) inhomogeneous distributions (see 
Kiyashchenko et al., 2003). The correlation dimension 
d(2) characterizes the degree of seismicity clustering. 
Lower values of d( 2) indicate stronger clustering in a 
wide range of scales. The value of amin characterizes the 
degree of seismicity clustering in seismically active 
parts of the considered region in some range of scales. 
The smaller values of cxmm can be interpreted as a 
manifestation of stronger clustering in densely 
populated areas. The value of ao is sensitive to the 
heterogeneity of the seismicity distribution in sparsely 
populated areas. The higher value of ao corresponds to 
more heterogeneous cases (see Kiyashchenko et al., 
2003). 


Results 

The results of calculation of multifractal characteristics 
for each particular EQ listed in Table 2 can be found in 
the paper by Kiyashchenko et al. (2003). Generally, 
the tendency of a decrease in the multifractal 
characteristics amin, S and d(2) and an increase in ao 
have been revealed before the mainshocks. As a case 
study, the case of well-known Hyogo-ken Nanbu 
(Kobe) EQ of 17 January, 1995 in Japan (cp = 34.583°N, 
A = 135.02°E, M = 7:2, depth =33 km) is presented in Fig. 
9. The length of sliding temporal window is taken 
equal to 0.5 x 10 8 s. It is seen that the decrease of the 
multifractal characteristics amin, S and d(2) and the 
increase of ao started approximately two years before 
the EQ. After the EQ, the multifractal characteristics 
tend to recover to their initial level. Such post-EQ 
behavior can be considered just as a recovery process, 
when the EQ focal system breaks down to the more 
disordered (chaotic) state after the release of main 
portion of seismic energy. 



G.OE+O l.QE+3 ZOE+8 3-OE+S 4.0E+8 
Time, s 

FIG. 9 DYNAMICS OF THE HIGHER-ORDER GENERALIZED 
FRACTAL DIMENSION (amin,), ENTROPY (S), CORRELATION 
DIMENSION d(2) AND ABSCISSA OF THE TOP OF 
MULTIFRACTAL SPECTRUM ao FOR A PARTICULAR CASE OF 
THE M=7.2 KOBE EQ WHICH OCCURRED IN JAPAN ON 
JANUARY 17, 1995. 

In order to highlight the general tendencies in the 
variations of multifractal characteristics prior to 
individual strong EQs, the superposed epoch method 
is applied to all the cases considered in Table 2. The 
procedure of the superposed epoch analysis is 
explained as follows. First, all the plots representing 
the variation of certain characteristics are shifted along 
the temporal axis, so that the moment t = 0 
corresponds to the main event occurrence time. Then, 
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the certain values of multifractal characteristics 
corresponding to a particular moment are stacked 
over all plots related to different EQs and averaged. 
The dynamics of the values ao, d(l)=S, d(2), amm, 
obtained by this method is shown in Fig. 10. 



Time, s 


FIG. 10 DYNAMICS OF THE SAME MULTIFRACTAL 
CHARACTERISTICS AS PRESENTED IN FIG. 9 FOR THE KOBE 
EQ, BUT FOR THE SUPERPOSED EPOCH SUMMARY OBTAINED 
FROM 12 EQS IN TABLE 2. 

For the superposed epoch analysis, the length of 
sliding temporal window for the calculation of 
multifractal characteristics is taken equal to 0.5 x 10 8 s 
for Japanese EQs and equal to 0.15 x 10 8 s for 
Californian EQs. The individual patterns of some EQs 
are neglected in this case, since the length of the 
temporal windows is not selected individually for 
each EQ. But even under such a condition, the general 
tendencies can be clearly seen in Fig. 10. So the 
multifractal characteristics OCmin , S and d(2), tend to 
decrease and ao tends to increase before the strong 
EQs. It can be treated as a signature of an increase in 
spatial inhomogeneity of the seismicity distribution in 
a wide range of scale levels, preceding strong EQs. 
Thus, it is found that the seismicity distribution 
evolves from the more homogeneous (disordered) 
state towards the more clustered (ordered) state. To 
compare the results for real and simulated seismicity, 
we have calculated the correlation dimension d(2), 
information dimension d(l), and the abscissa of the top 
of multifractal spectrum ao in a sliding temporal 
window of 20000 s for the synthetic seismicity catalog, 
obtained as the result of simulation of the destruction 
process (see Kiyashchenko et al. (2004) for details). The 
results are shown in Fig. 11. The calculations of the 
high-order generalized fractal dimension amm, have not 
been performed because of too small number of events. 
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It is clearly seen that the correlation and information 
dimensions tend to decrease before the main rupture 
moment. The multifractal characteristics ao increase 
before the main rupture. These results support the 
conclusion concerning the increase of heterogeneity of 
the seismicity distribution in a wide range of scales 
prior to main rupture. So the tendencies which have 
been revealed in the behavior of real seismicity are in 
agreement with the results of simulation of the 
destruction of elastic body with a number of shear 
cracks reported in Kiyashchenko et al. (2004). 
Therefore, it can be argued that multifractal analysis of 
the seismicity distribution dynamics is one of possible 
sources of information on the precursory evolution of 
the Earth's crust. 



FIG. 11 DYNAMICS OF MULTIFRACTAL CHARACTERISTICS OF 
THE SYNTHETIC SEISMICITY DISTRIBUTION OBTAINED AS A 
RESULT OF SIMULATION OF DESTRUCTION PROCESS. 

Discussion and conclusions 

We have analyzed the world-wide studies of 
multifractal properties of seismicity and highlighted 
certain tendencies in dynamics of the seismicity 
multifractal spectra prior to main strong EQs. As it 
follows from the corresponding researches (Goltz, 
1997), the seismicity multifractal spectrum extends in 
such a way that amin and a max, shift to the right and ao 
increases prior to strong EQs. Those peculiarities mean 
the transition from homogeneous (random) to 
heterogeneous (ordered, clustered) patterns; 
intensification of the clustering within rupture areas; 
increment of the local fractal dimension; and a change 
in the ratio of highly clustered and sparsely populated 
areas. However, those results have been obtained on 
the basis of poor statistics. So we have fulfilled an 
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additional multifractal analysis of the seismicity 
distributions before 12 strong EQs occurring in Japan 
and Southern California during an extended period of 
more than 10 years (see Table 2). We have tried to 
receive sufficient statistics to draw a definite 
conclusion about seismicity multifractal dynamics, 
and have calculated four multifractal characteristics of 
seismicity, namely, the minimal value of the Holder 
exponent OCmin — d(+oo), the information dimension S = 
ai = d(l), the correlation dimension d( 2), and the 
abscissa of the top of multifractal spectrum ao, for 
seismicity distribution around 12 strong EQs listed in 
Table 2. Then their dynamics in the process of 
approach of the main EQ has been studied. A clear 
tendency of decrement of the multifractal 
characteristics OCmin, S and d(2), and increment of ao 
before the main shocks have been revealed. It can be 
interpreted as a signature of increase of spatial 
inhomogeneity of the seismicity distribution in a wide 
range of scale levels, preceding strong EQs. Thus, 
seismicity distribution evolves from the more 
homogeneous (disordered) state towards the more 
clustered (ordered) state. That is consistent with the 
few world-wide studies. Also the revealed dynamics is 
consistent with the dynamics of multifractal 
characteristics of the synthetic seismicity distribution 
obtained as a result of the destruction process 
simulation (Kiyashchenko et al., 2004). In view of our 
reliable statistics, it is possible to argue that the 
variations of multifractal characteristics presented here 
appear real due to the processes of the reorganization 
of the seismicity but not due to numerical errors or 
instability of the procedure of calculation of 
generalized fractal dimensions. So, the results 
obtained offer an opportunity to use multifractal 
scaling characteristics of the seismicity for monitoring 
of the destruction process in seismically active areas 
and thus for a forecast of strong seismic events. 
Appreciable information on pre-rupture evolution of 
fault (crack) network in the lithosphere could be 
provided by the study of such multifractal dynamics. 

On the basis of the research fulfilled in this paper, it 
can be concluded that multifractal analysis of 
geophysical data could be a promising tool for 
extraction of the precursory signatures of the extreme 
natural events including destructive EQs. However, 
the anomalous variations in multifractal parameters of 
seismicity distribution is found to begin a few years 
before the main shock, which seems to be consistent 
with the former corresponding results by Haikun et al. 
(1993) and Hirabayashi et al. (1992). Thus, this kind of 
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EQ precursor based on the seismicity distribution is 
likely to be a typical medium-term EQ prediction. 
When any particular region with expectation of a big 
EQ has been known from the seismological point of 
view, the method proposed in this paper would be of 
extreme importance in predicting EQ. In order to 
realize the real EQ prediction, it is highly required to 
have additionally the corresponding analysis for 
electromagnetic effects (Molchanov and Hayakawa, 
2008; Hayakawa and Ida, 2008)), because these 
electromagnetic effects are considered to be promising 
candidates of short-term precursors. 
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